/*
||
|| adv.U
||
|| This file contains routines that are used to take the
|| solution through one time step.  Some of the required 
|| reservoir code is in resadv.u.  Adapted from a CRAY 2
|| code, this code has been radically modified to take
|| advantage of the ANL macro set programming model.
|| Reservoir update is also handled in this module.
||
*/

/* Macro for ANL declarations */
EXTERN_ENV


/* include all mp3d data defintions */
#define MAIN extern
#include <stdio.h>
#include "common.h"

/*
|| BOUNDARY :  Takes passed molecules and performs the
|| required modifications to its position and velocity vector based
|| on the boundary condition applicable to the current cell.
|| Boundry conditions at a particular cell are represented as
|| equations of a plane in that cell in a surface normal
|| representation.  Boundry temperature and absorbtion coeffecient
|| are also stored for handling diffuse (non-specular) boundary
|| conditions.
*/

static int 
  local_BC[MAX_PROCESSORS],              /* boundary collision counter */
  local_num_collisions[MAX_PROCESSORS];  /* molecule collision counter */

static void boundary(int bn, struct particle *Part, int I)
{
  int local_next;
  float d2,
        temp;
  struct bcond *BCond = &Bconds[bn];
  struct ares *RPart;

  ++local_BC[I];
  if (bn < COND_OFFSET)  /* Solid wall case */ {
    d2 = 2*(BCond->bc_A*Part->x 
          + BCond->bc_B*Part->y 
          + BCond->bc_C*Part->z 
          + BCond->bc_D);
    if (d2 < 0.0) {
      Part->x -= d2 * BCond->bc_A;
      Part->y -= d2 * BCond->bc_B;
      Part->z -= d2 * BCond->bc_C;
      Part->where = &Cells[CELL(Part->y)][CELL(Part->z)][CELL(Part->x)];
      d2 = 2*(BCond->bc_A * Part->u 
            + BCond->bc_B * Part->v
            + BCond->bc_C * Part->w);
      temp = BCond->bc_temp*(1.0 - BCond->bc_coef);
      Part->u = (Part->u - d2*BCond->bc_A)
        * BCond->bc_coef + BCond->bc_A*temp;
      Part->v = (Part->v - d2*BCond->bc_B)
        * BCond->bc_coef + BCond->bc_B*temp;
      Part->w = (Part->w - d2*BCond->bc_C)
        * BCond->bc_coef + BCond->bc_C*temp;
    }

  } else /* entrance or exit condition */ {
    bn -= COND_OFFSET;
    switch (bn) {
      case 1:  /* exit */
        LOCK(Global->next_res_lock)
          local_next = Global->next_res++ % num_res;
          Global->exit_num_flow++;
        UNLOCK(Global->next_res_lock)
        RPart = &Ares[local_next];
        Part->x = 2.0 + RPart->xr*fluid_speed;
        Part->y = 2.0 + RPart->yr*space_heightm4;
        Part->z = 2.0 + RPart->zr*space_widthm4;
        Part->where = &Cells[CELL(Part->y)][CELL(Part->z)][CELL(Part->x)];
        Part->u = RPart->ur + fluid_speed;
        Part->v = RPart->vr;
        Part->w = RPart->wr;
        Part->r = RPart->rr;
        Part->s = RPart->sr;
        break;

      case 2:  /* entrance */
        LOCK(Global->next_res_lock)
          local_next = Global->next_res++ % num_res;
        UNLOCK(Global->next_res_lock)
        RPart = &Ares[local_next];
        Part->x = 2.0 + RPart->xr*fluid_speed;
        Part->y = 2.0 + RPart->yr*space_heightm4;
        Part->z = 2.0 + RPart->zr*space_widthm4;
        Part->where = &Cells[CELL(Part->y)][CELL(Part->z)][CELL(Part->x)];
        Part->u = RPart->ur + fluid_speed;
        Part->v = RPart->vr;
        Part->w = RPart->wr;
        Part->r = RPart->rr;
        Part->s = RPart->sr;
        break;

      default:
        LOCK(Global->output)
          printf("ERROR: Unknown Condition = %d.\n", bn);
        UNLOCK(Global->output)
        MAIN_END(-1)
        exit(1);
    } /* switch */
  } /* if condition */
} /* boundary() */

/*
|| ADD TO FREE STREAM :  Adds molecules to entrance to build
|| up to the already calculated free stream density.
*/
static void add_to_free_stream (int I) 
{
  int i,             /* subcript for dist loop */
      local_next,
      local_num_mol;
  struct ares *RPart;
  struct particle *Part;

  /* Must synchronize here to let master calc num to add */
  BARRIER(Global->end_sync, number_of_processors)

  /* Add entrance molecules up to free flow number */
  /* as calculated in master advance               */
  for (i = I; i < Global->max_sub; i+= number_of_processors) { 
    LOCK(Global->num_mol_lock)
      local_num_mol = Global->num_mol++;
    UNLOCK(Global->num_mol_lock)

    if (local_num_mol < MAX_MOL) {
      LOCK(Global->next_res_lock)
        local_next = Global->next_res++ % num_res;
      UNLOCK(Global->next_res_lock)

      Part = &Particles[local_num_mol];
      RPart = &Ares[local_next];
      Part->x = 2.0 + RPart->xr * fluid_speed;
      Part->y = 2.0 + RPart->yr * space_heightm4;
      Part->z = 2.0 + RPart->zr * space_widthm4;
      Part->where = &Cells[CELL(Part->y)][CELL(Part->z)][CELL(Part->x)];
      Part->u = RPart->ur + fluid_speed;
      Part->v = RPart->vr;
      Part->w = RPart->wr;
      Part->r = RPart->rr;
      Part->s = RPart->sr;

    } else {
      LOCK(Global->output)
        printf("ERROR: Out of molecule space.\n");
      UNLOCK(Global->output)
      MAIN_END(-1)
      exit(1);
    }
  }
  BARRIER(Global->end_sync, number_of_processors)
} /* add_to_free_stream() */

/*
|| RESET :  This routine prepares the space for the next time-step. 
|| Only positive numbers are zeroed in the space, as negative ones
|| represent boundarys.
*/
static void reset (int I) 
{
  int i,
      j,
      k,
      iter,
      loop_max,
      nz;
  struct acell *ACell;

  /* reset space */
  nz = num_zcell - 4;
  loop_max = (num_ycell - 4)*nz;
  for (iter = I; iter < loop_max; iter += number_of_processors) {
    j = 2 + iter/nz;
    k = 2 + iter%nz;
    for (i = 0; i < num_xcell; i++) {
      ACell = &Cells[j][k][i];
      ACell->avg_prob = 0.99*ACell->avg_prob 
                      + 0.01*Global->coll_prob
                      * (float) ACell->cell_pop;
      ACell->cell_pop = 0;
      if (ACell->space > 0)
        ACell->space = 0;
    }
  } 

  /* Reset space lying partially under boundary conditions */
  for (i = I; i < bc_count; i+= number_of_processors) {
    Bc_space[i] = 0;
  }
  BARRIER(Global->end_sync, number_of_processors)
} /* reset() */

/*
|| COLLIDE:  Collides the passed molecules in a collision of
|| statistically random outcome.
*/
static void collide(struct particle *Part1, struct particle *Part2, int I)
{
/* I is processor id */
  int *pm;    /* pointer into mix table */
  float *ps,  /* pointer into sign table */
        r_mean,
        rel[5],
        s_mean,
        temp,
        u_mean,
        v_mean,
        w_mean;

  /* Calculate relative and mean velocities */
  rel[0]  = (Part1->u - Part2->u)*0.5;
  u_mean = Part1->u - rel[0];

  rel[1]  = (Part1->v - Part2->v)*0.5;
  v_mean = Part1->v - rel[1];

  rel[2]  = (Part1->w - Part2->w)*0.5;
  w_mean = Part1->w - rel[2];

  rel[3]  = (Part1->r - Part2->r)*0.5;
  r_mean = Part1->r - rel[3];

  rel[4]  = (Part1->s - Part2->s)*0.5;
  s_mean = Part1->s - rel[4];

  /* Lookup random collision outcome */
  pm = Mix_table[TCC_random() % MIXS_POSSIBLE];
  ps = Sign_table[TCC_random() % SIGNS_POSSIBLE];

  /* perform collision */
  temp = (*ps) * rel[*pm];
  Part1->u = u_mean + temp;
  Part2->u = u_mean - temp;
  ps++;
  pm++;
  temp = (*ps) * rel[*pm];
  Part1->v = v_mean + temp;
  Part2->v = v_mean - temp;
  ps++;
  pm++;
  temp = (*ps) * rel[*pm];
  Part1->w = w_mean + temp;
  Part2->w = w_mean - temp;
  ps++;
  pm++;
  temp = (*ps) * rel[*pm];
  Part1->r = r_mean + temp;
  Part2->r = r_mean - temp;
  ps++;
  pm++;
  temp = (*ps) * rel[*pm];
  Part1->s = s_mean + temp;
  Part2->s = s_mean - temp;

  ++local_num_collisions[I];
} /* collide() */

static void move_single(int pn, int I)
{
  float u,
        v,
        w;
  int bc_val,
      space_value;
  struct acell *ACell;
  struct particle *Part;

  Part = &Particles[pn];
  ACell = Part->where;

  /* Do sample from previous time-step */
  u = Part->u;
  v = Part->v;
  w = Part->w;
#ifdef LOCKING
  LOCK(ACell->sample_lock)
#endif
    ACell->U += u;
    ACell->V += v;
    ACell->W += w;
    ACell->Usqr += u*u;
    ACell->Vsqr += v*v;
    ACell->Wsqr += w*w;
#ifdef LOCKING
  UNLOCK(ACell->sample_lock)
#endif

  /* Move them x=x+u */
  Part->x += u;
  Part->y += v;
  Part->z += w;
  ACell = &Cells[CELL(Part->y)][CELL(Part->z)][CELL(Part->x)];
  Part->where = ACell;

  /* boundary interaction if required */
#ifdef LOCKING
  LOCK(ACell->space_lock)
#endif
    space_value = ACell->space;
    if (space_value < 0) {
#ifdef LOCKING
    UNLOCK(ACell->space_lock)
#endif
    boundary(-space_value, Part, I);
    ACell = Part->where;
#ifdef LOCKING
    LOCK(ACell->space_lock)
#endif
    }

    /* Update cell populations */
    ACell->cell_pop++;
    ACell->cum_cell_pop++;

    /* Do collision pairing */
    space_value = ACell->space;
    if (space_value == 0) {
      ACell->space = pn;

    } else {
      if (space_value > 0) {
        ACell->space = 0;
        if (RANDOM <= ACell->avg_prob)
          collide(&Particles[space_value], Part, I);

      } else { /* must use boundary space */
        if (-space_value < COND_OFFSET) {
          bc_val = Bc_space[-space_value];
          if (bc_val != 0) {
            if (RANDOM <= ACell->avg_prob)
              collide(&Particles[bc_val], Part, I);
          } else
            Bc_space[-space_value] = pn;
        }
      }
    }
#ifdef LOCKING
  UNLOCK(ACell->space_lock)
#endif
} /* move_single() */

/*
|| MOVE :  Move Molecules through one time-step, doing
|| boundary interactions collisions as required.  A processor
|| be it master or slave is assigned to work on one molecule
|| at a time via a distributed loop.
*/
static void move (int I) 
{
  int j,
      limit,
      pn,
      step;

  local_num_collisions[I] = 0;
  local_BC[I] = 0;

  step = number_of_processors*clump;
  limit = clump*(Global->num_mol/clump);
  for (j = I*clump; j < limit; j += step)
    for (pn = j; pn < j + clump; pn++)
      move_single(pn,I);

  if (I == number_of_processors - 1)
    for (pn = limit; pn < Global->num_mol; pn++) 
      move_single(pn,I);

  LOCK(Global->collisions_lock)
    Global->num_collisions += local_num_collisions[I];
    Global->total_collisions += local_num_collisions[I];
    Global->BC += local_BC[I];
  UNLOCK(Global->collisions_lock)
  BARRIER(Global->end_sync, number_of_processors)
} /* move() */

/*
|| RESERVOIR MOVE :  Moves mollecules in the reservoir through 
|| a single time-step and performs cyclic boundary conditions.
*/
static void reservoir_move (int I) 
{
  int pn;

  struct ares *RPart;

  for (pn = I; pn < num_res; pn += number_of_processors) {

    /* Move them and keep in unit cube by employing 
       cyclic boundary conditions */
    RPart = &Ares[pn];
    RPart->xr += RPart->ur + 2.0;
    RPart->xr -= (int) RPart->xr;

    RPart->yr += RPart->vr + 2.0;
    RPart->yr -= (int) RPart->yr;

    RPart->zr += RPart->wr + 2.0;
    RPart->zr -= (int) RPart->zr;
  }
  BARRIER(Global->end_sync, number_of_processors)
} /* reservoir_move() */

/*
|| RESERVOIR COLLIDE:  Performs collision of each reservoir molecule
|| on every time-step.  Uses different pairings for each
|| time-step to keep the whole system mixed.
*/
static void reservoir_collide (int I) 
{
  float rel[5],
        u_mean,
        v_mean,
        w_mean,
        r_mean,
        s_mean,
        *ps,  /* pointer into sign table */
        temp;
  int *pm,  /* pointer into mix table */
      pn,
      pn1,
      pn2;
  struct ares *RPart1,
              *RPart2;

  for (pn = I; pn < num_res/2; pn += number_of_processors) {
    pn1 = 2*pn; /* even */
    pn2 = (pn1 + Global->res_offset) % num_res; /* odd */
    RPart1 = &Ares[pn1];
    RPart2 = &Ares[pn2];

    /* Get mean and relative velocity */
    rel[0]  = (RPart1->ur - RPart2->ur)*0.5;
    u_mean = RPart1->ur - rel[0];

    rel[1]  = (RPart1->vr - RPart2->vr)*0.5;
    v_mean = RPart1->vr - rel[1];

    rel[2]  = (RPart1->wr - RPart2->wr)*0.5;
    w_mean = RPart1->wr - rel[2];

    rel[3]  = (RPart1->rr - RPart2->rr)*0.5;
    r_mean = RPart1->rr - rel[3];

    rel[4]  = (RPart1->sr - RPart2->sr)*0.5;
    s_mean = RPart1->sr - rel[4];

    /* Pick a random collision outcome */
    pm = Mix_table[TCC_random() % MIXS_POSSIBLE];
    ps = Sign_table[TCC_random() % SIGNS_POSSIBLE];

    /* Do the collision */
    RPart1->ur = u_mean + (temp = (*ps) * rel[*pm]);
    RPart2->ur = u_mean - temp;
    ps++;
    pm++;
    RPart1->vr = v_mean + (temp = (*ps) * rel[*pm]);
    RPart2->vr = v_mean - temp;
    ps++;
    pm++;
    RPart1->wr = w_mean + (temp = (*ps) * rel[*pm]);
    RPart2->wr = w_mean - temp;
    ps++;
    pm++;
    RPart1->rr = r_mean + (temp = (*ps) * rel[*pm]);
    RPart2->rr = r_mean - temp;
    ps++;
    pm++;
    RPart1->sr = s_mean + (temp = (*ps) * rel[*pm]);
    RPart2->sr = s_mean - temp;
  }
  BARRIER(Global->end_sync, number_of_processors)
} /* reservoir_collide() */
  
/*
|| SLAVE ADVANCE SOLUTION :  Advances the solution take_steps 
|| time-steps for the molecules All Boundry conditions and 
|| collisions are performed.  This routine represents the code
|| executed by all slave processors (all but one).  MASTER ADVANCE
|| SOLUTION is the same execpt it also increments take_steps after
|| each steps.
*/
void slave_advance_solution() {
  int I;

  GET_PID(I)
  while (Global->take_steps>0) {
    reset (I);              /* reset the memory map of space */
    move (I);               /* move molecules one delta t and collide */
    add_to_free_stream (I); /* add to free stream */
    reservoir_move (I);     /* move reservoir molecules */
    reservoir_collide (I);  /* collide reservoir molecules */
  }
} /* slave_advance_solution() */

/*
|| MASTER ADVANCE SOLUTION :  Advances the solution take_steps 
|| time-steps for the molecules All Boundry conditions and 
|| collisions are performed.  This routine represents the code
|| executed by the host  processor.  MASTER ADVANCE SOLUTION is 
|| the same as SLAVE ADVANCE SOLUTION except it also increments 
|| take_steps after each steps.  Remember explicit barriers are not
|| required as distributed loops in other routines imposed
|| barriers.  Master advance also does small amounts of work that
|| must not be duplicated by slaves.
*/
void master_advance_solution() {
  int I;

  GET_PID(I)
  while (Global->take_steps>0) {
    /* Reset counters. */
    Global->BC = 0;
    Global->num_collisions = 0;

    reset (I);              /* reset map of space */
    Global->take_steps--; /* safe to decrement number of steps */

    if (Global->take_steps == 2) {
#ifdef AUG
       ST_LOG(ST_RESET_STATS);
#else
       RESET_STATS
#endif
    }
    move (I);               /* move through delta t */

    /* How many new molecules should we add? */
    Global->tot_flt_input += Global->num_flow;
    Global->tot_input = (long) (Global->tot_flt_input);
    Global->max_sub = Global->tot_input - Global->exit_num_flow;
    Global->res_offset += 2;

    add_to_free_stream (I); /* add new molecules */

    /* Adjust entrance counts */
    if (Global->exit_num_flow < Global->tot_input)
      Global->exit_num_flow = Global->tot_input;

    reservoir_move (I);     /* move reservoir molecules */
    reservoir_collide (I);  /* collide reservoir molecules */
  } /* while more steps */
} /* master_advance_solution() */
